function diff = q_1_fun(m,param)

    if param.zeta==0
        diff = ones(size(m));
    else
        f = @(x) ( (1-param.alpha).*2./param.sigma^2.*param.zeta./x./q_0_fun(x,param));
        diff = nan(size(m));
        for i=1:length(m)
            diff(i)= param.psi - integral(f,param.m_l,m(i));
        end
    end

end